EstuarineFishGillMicrobiota_part2_21.07.24
2024-07-21
here we perform some overall alpha diversity analyses inspecting species richness and diversity as well as determining the fish gill mucus core microbiota via prevalence and abundance. Spatio-temporal of bacterioplankton and fish is found in chapter 4.5 depicting especially the interesting microbiota adapation in migrating fish. First steps of data exploration for modeling alpha diversity are shown in 5. however seem to transfer low information.
1 Global Setup
1.1 Path
before you start: .rs.restartR()
1.2 Packages
1.3 Functions
1.4 Input Files
1.5 Tutorials
1.6 Setup Analysis
#-
4 Alpha Overview
4.1 Observed Diversity Plots
simple richness plot over all samples
##########################
#Observed Diversity Plots#
##########################
for (i in names(pslist)[grepl("ps", names(pslist))]){
require(plyr); require(ggrepel); require(cowplot); require(phyloseq)
#if (OperatingSystem == "Mac" ) {quartz() }
tryCatch({
g <- paste(i); print(i)
#gg <- sub('ps', '', g)
A<- plot_richness(pslist[[i]], x = "SampleID", measures = c("Observed")) + # Shannon
geom_bar(aes(fill = sample_data(pslist[[i]])$Species), stat = "identity", position = "stack") + #fill = Phylum
scale_x_discrete(limits = SSU_Samples[SSU_Samples %in% sample_names(pslist[[i]])]) +
#scale_x_discrete(breaks=factor(sample_data(pslist[[i]])$Replicates, levels=SLOrder)) +
scale_fill_manual(values = col.Palette$col.Palette.Species) +
theme(axis.text.x.bottom = element_text(size=rel(0.8),angle = 45, vjust = 1, hjust = 1),
legend.title = element_text( size = 6),
legend.text = element_text(size = 6),
legend.key.size = unit(0.4, 'cm'),)+
labs(x = element_blank(), y = "Alpha Diversity Measure (non normalized)")
title <- ggdraw() + draw_label_themeRK(paste(g), element = "plot.title",x = 0.05, hjust = 0, vjust = 1)
subtitle <- ggdraw() + draw_label_themeRK(paste("Paired End 16S",sep = " "), element = "plot.subtitle",
x = 0.05, hjust = 0,
vjust = 1)
prow <- cowplot::plot_grid(A, labels = c(""), ncol = 1)
A<- plot_grid(title, subtitle, prow, ncol = 1, rel_heights = c(0.04, 0.05, 0.98))
plot(A)
ggsave(A, filename = paste(paste(save_name, Type, paste(g, "Observed-AlphaDiversity", sep="_"), Date, sep="_"), ".png", sep=""), path = pathPlots, device='png', dpi=300, width = 8,
height = 8)
}, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})}## [1] "ps_SSU"
## [1] "ps_WF"
## [1] "ps_Fish"
## [1] "ps_OE"
## [1] "ps_GC"
## [1] "ps_OEWF"
## [1] "ps_GCWF"
4.2 Species Accumulation
Visualizing of diversity is covered in raw and filtered datasets
for (Dataset in names(pslistraw)) {
require(vegan)
spacc_raw_ps <- pslistraw[[Dataset]]
spacc_ps <- pslist[[Dataset]]
Datasetname <- sub("ps_", "", Dataset)
spacc_raw <- vegan::specaccum(otu_table(spacc_raw_ps), method = "random", permutations = 100)
spacc <- vegan::specaccum(otu_table(spacc_ps), method = "random", permutations = 100)
#creating a dataframe for ggplot2
data_spacc_raw <- data.frame(Specimen=spacc_raw$sites, Richness=spacc_raw$richness, SD=spacc_raw$sd)
A<- ggplot() +
geom_point(data=data_spacc_raw, aes(x=Specimen, y=Richness)) +
geom_line(data=data_spacc_raw, aes(x=Specimen, y=Richness)) +
geom_ribbon(data=data_spacc_raw ,aes(x=Specimen, ymin=(Richness-2*SD),ymax=(Richness+2*SD)),alpha=0.2) +
atheme +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
#legend.background = element_rect(fill='transparent'), #transparent legend bg
#legend.box.background = element_rect(fill='transparent') #transparent legend panel
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=10),
legend.text=element_text(size=10)) +
theme(
panel.grid.major = element_line(colour = "grey50"),
panel.grid.minor = element_line(colour = "grey50")) +
ggtitle(paste(Datasetname, "raw"))
Species_Accumulation_raw <- cowplot::plot_grid(A, labels = c(""), ncol = 1)
data_spacc <- data.frame(Specimen=spacc$sites, Richness=spacc$richness, SD=spacc$sd)
B<- ggplot() +
geom_point(data=data_spacc, aes(x=Specimen, y=Richness)) +
geom_line(data=data_spacc, aes(x=Specimen, y=Richness)) +
geom_ribbon(data=data_spacc,aes(x=Specimen, ymin=(Richness-2*SD),ymax=(Richness+2*SD)),alpha=0.2) +
atheme +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
#legend.background = element_rect(fill='transparent'), #transparent legend bg
#legend.box.background = element_rect(fill='transparent') #transparent legend panel
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=10),
legend.text=element_text(size=10)) +
theme(
panel.grid.major = element_line(colour = "grey50"),
panel.grid.minor = element_line(colour = "grey50")) +
ggtitle(paste(Datasetname, "abundance filtered"))
Species_Accumulation_Filtered <- cowplot::plot_grid(B, labels = c(""), ncol = 1)
require(vegan)
cowplot::plot_grid(Species_Accumulation_raw, Species_Accumulation_Filtered, labels = c("A", "B"), ncol = 1) -> part_1
ggsave(part_1, filename = paste(paste(save_name, Type, "SpeciesAccumulation", Datasetname, Date, sep="_"), ".png", sep=""), path = pathPlots , device='png', dpi=300, width = 6,height = 9)
plot(part_1)
}
4.3 Alpha Diversity Barplot
simple plotting of the taxa with highest abundance over time and space clue here are the preserved and automatically picked colors between taxa
##########################
#Barplot with color code#
##########################
#Run once to determine color for both fish species
TaxLevel <- "Genus"
for (Dataset in names(pslist)[grepl("ps_Fish", names(pslist))]){
require(plyr); require(dplyr); require(ggrepel); require(cowplot); require(phyloseq)
Datasetname <- sub("ps_", "", Dataset)
FILENAME <- paste(paste(save_name, "Alpha_BarPlot_Publication", TaxLevel, Datasetname, sep="_"), ".png", sep="")
WIDTH = 10 + length(sample_names(pslist[[Dataset]])) *0.12
if (Datasetname %in% c("GC", "GCWF", "Fish")) {
pslist[[Dataset]] <- prune_samples(sample_names(pslist[[Dataset]])[!grepl("GCSU21BBEB6", sample_names(pslist[[Dataset]]))], pslist[[Dataset]])
} else {pslist[[Dataset]] <- pslist[[Dataset]]}
################################
#Create Relative Abundance Data#
################################
ps_alpha_barplot <- pslist[[Dataset]] %>%
tax_glom(taxrank = TaxLevel) %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>% # Transform to rel. abundance
psmelt() %>% # Melt to long format
filter(Abundance > 1) %>% # Filter out low abundance taxa
arrange(Genus) %>% # Sort data frame alphabetically by phylum
dplyr::arrange(desc(Abundance))
Order <- SSU_Samples[SSU_Samples %in% sample_names(pslist[[Dataset ]])]
ps_alpha_barplot$ASV <- sub(".*:", "", ps_alpha_barplot$OTU)
ps_alpha_barplot$ASV <- sub('\\.', ' ', ps_alpha_barplot$ASV)
############################
#Create TotalAbundance Data#
############################
phylum_abundance <- ps_alpha_barplot %>%
dplyr::group_by(.data[[TaxLevel]]) %>%
dplyr::summarise(TotalAbundance = sum(Abundance))
ordered_levels <- phylum_abundance %>%
dplyr::arrange(desc(TotalAbundance)) %>%
pull(.data[[TaxLevel]])
ps_alpha_barplot$Taxa <- factor(ps_alpha_barplot[[TaxLevel]], levels = ordered_levels)
ps_alpha_barplot$SampleID <- factor(ps_alpha_barplot$SampleID, levels = SSU_Samples)
##########################
#Subset for Interest ASVs#
##########################
df <- ps_alpha_barplot #[ps_alpha_barplot$OTU %in% GroupOfInterest,]
######################
#Define Phylum Colors#
######################
phylum_colors <- c(
"Actinobacteriota" = "red",
"Cyanobacteria" = "green",
"Bacteroidota" = "darkgreen",
"Proteobacteria1" = "darkorange2",
"Proteobacteria2" = "#663300",
#"Patescibacteria" = "pink",
"Acidobacteriota " = "#FF00CC",
"Planctomycetota" = "purple",
"Verrucomicrobiota" = "blue",
"Chloroflexota" = "cyan",
"Deinococcota" = "#FFFF00",
"Gemmatimonadota" = "#8F7C00",
"Other" = "grey20")
generate_shades <- function(base_color, num_variations) {
color_ramp <- colorRampPalette(c(base_color, "white"))
# Generate a vector of colors
colors <- color_ramp(num_variations+1)
# Create a data frame with the colors
color_palette <- colors[1:(length(colors) - 1)]
return(color_palette)}
##################################################
#Create Color Dataframe for Taxa from PhylaColors#
##################################################
taxonomy_df <- df[c("Phylum", "Taxa")]
#taxonomy_df$ASV <- sub(".*:", "", taxonomy_df$OTU)
taxonomy_df <- taxonomy_df[!duplicated(taxonomy_df$Taxa),]
#############################################
#Order to ordered_levels from Total Abudance#
order_index <- match(ordered_levels, taxonomy_df$Taxa)
taxonomy_df <- taxonomy_df[order_index, ]
print(head(taxonomy_df[c("Phylum", "Taxa")]))
taxonomy_df$Phylum2 <- ifelse(taxonomy_df$Phylum == "Proteobacteria",
paste0("Proteobacteria", ave(seq_along(taxonomy_df$Phylum),
taxonomy_df$Phylum, FUN = function(x) ifelse(seq_along(x) %% 2 == 0, 2, 1))),
taxonomy_df$Phylum)
unique_phyla <- unique(taxonomy_df$Phylum2)
# Assign base colors to unique phyla
#phylum_colors <- setNames(base_colors[1:length(unique_phyla)], unique_phyla)
# Print the mapping of phyla to colors
phylum_cols<- as.data.frame(phylum_colors)
phylum_cols$Phyla <- rownames(phylum_cols)
colored_taxonomy <- data.frame(
Taxa = taxonomy_df$Taxa,
Phylum2 = taxonomy_df$Phylum2,
Color = character(length(taxonomy_df$Taxa)), # Initialize an empty character vector
stringsAsFactors = FALSE)
# Generate shades and assign colors to each taxon
for (ii in seq_along(unique_phyla)) {
phylum <- unique_phyla[ii]
if (phylum %in% names(phylum_colors)) {
base_color <- phylum_colors[[phylum]]
} else {
base_color <- phylum_colors[["Other"]]
}
taxon_indices <- taxonomy_df$Phylum2 == phylum
num_variations <- sum(taxon_indices)
shades <- generate_shades(base_color, pmin(num_variations, 5))
#colored_taxonomy$Color[taxon_indices] <- shades
colored_taxonomy$Color[which(taxon_indices)[1:pmin(num_variations, 5)]]<- shades
}
colored_taxonomy$Color <- gsub("^$", NA, trimws(colored_taxonomy$Color))
colored_taxonomy$Color <- replace_na(colored_taxonomy$Color, "grey50")
##################
# #Some adjustments#
updateColor <- function(taxa, color) {
if (any(colored_taxonomy$Taxa == taxa)) {
colored_taxonomy[colored_taxonomy$Taxa == taxa,]$Color <- color
}
return(colored_taxonomy)
}
colored_taxonomy <- updateColor("Psychrobacter", "#FFCC00")
colored_taxonomy <- updateColor("Polynucleobacter", "#FF6600")
#updateColor("Photobacterium", "darkorange3")
# updateColor("Vibrio", "#FF9966")
#updateColor("Citrobacter", "#FFCC00")
#updateColor("Psychrobacter", "#FFCC00")
# updateColor("Enterobacter", "#993300")
#updateColor("Acinetobacter", "#FF9933")
Alpha_Diversity_df <- left_join(df, colored_taxonomy)
color_mapping <- setNames(Alpha_Diversity_df$Color, Alpha_Diversity_df$Taxa)
levels(Alpha_Diversity_df$Taxa)[!levels(Alpha_Diversity_df$Taxa) %in% ordered_levels[1:34]] <-
"Other"
p <- ggplot(Alpha_Diversity_df,
aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
facet_grid(. ~ factor(Reps, levels = RepOrder), drop = TRUE, scale = "free", space = "free_x") +
atheme2 +
ylab("Relative Abundance [%] \n") + xlab("") +
labs(fill="") +
scale_fill_manual(values = color_mapping) +
#ggrepel::geom_text_repel(x = df$SampleID, y = df$Abundance, aes(label = Labels),
#inherit.aes = FALSE, min.segment.length = 0,nudge_y = 20,size=3, max.overlaps = Inf) +
#awhite +
theme(strip.text.y = element_text(angle = 0)) +
#theme(axis.text.x=element_blank(),
#axis.text.y=element_blank(),
#axis.title.y.left = element_blank(),
#axis.ticks.y =element_blank()
#) +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
#legend.background = element_rect(fill='transparent'), #transparent legend bg
#legend.box.background = element_rect(fill='transparent') #transparent legend panel
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=9),
legend.text=element_text(size=9),
axis.text.x.bottom = element_text(size= 9, face = "bold", angle = 45, hjust = 1),
strip.text.y.left = element_text(size = 9,face = "bold"),
axis.title.y.left = element_text(size=12,face = "bold"))
p
COL <- col.Palette$col.Palette.Reps[names(col.Palette$col.Palette.Reps) %in% RepOrder]
g <- ggplot_gtable(ggplot_build(p))
stripr <- which(grepl('strip-t', g$layout$name))
fills <- alpha(COL, 0.8)
k <- 1
for (iii in stripr) {
j <- which(grepl('rect', g$grobs[[iii]]$grobs[[1]]$childrenOrder))
g$grobs[[iii]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1}
Alpha_Diversity_BarPlot <- cowplot::plot_grid(g, labels = c(""), ncol = 1)
Alpha_Diversity_BarPlot <- plot_grid(Alpha_Diversity_BarPlot, ncol = 1, rel_heights = c(100))
ggsave(Alpha_Diversity_BarPlot, filename = FILENAME, path = pathPlots, device='png', dpi=300, width = WIDTH, height = 6)
}## Phylum Taxa
## 10 Bacteroidota Elizabethkingia
## 204 Proteobacteria Enterobacter
## 196 Proteobacteria Citrobacter
## 12 Proteobacteria Psychrobacter
## 218 Bacteroidota Chryseobacterium
## 53 Verrucomicrobiota Luteolibactercolored_taxonomy_Fish <- colored_taxonomy
pslist[["Optics"]] <- list()
pslist[["Optics"]][["colored_taxonomy_Fish"]] <- colored_taxonomy_Fish
##############################
#Now use color from both fish#
##############################
TaxLevel <- "Genus"
for (Dataset in names(pslist)[grepl("ps_OE|ps_GC|ps_WF", names(pslist))]){
require(plyr); require(dplyr); require(ggrepel); require(cowplot); require(phyloseq)
Datasetname <- sub("ps_", "", Dataset)
FILENAME <- paste(paste(save_name, "Alpha_BarPlot_Publication", TaxLevel, Datasetname, sep="_"), ".png", sep="")
if (Datasetname %in% c("GC", "OE")) {
WIDTH <- 10 + length(sample_names(pslist[[Dataset]])) *0.12
} else if (Datasetname %in% c("WF")) {
WIDTH <- 5 + length(sample_names(pslist[[Dataset]])) *0.12
} else {
WIDTH <- 10 + length(sample_names(pslist[[Dataset]])) *0.12
}
ps <- pslist[[Dataset]]
if (Datasetname %in% c("GC", "GCWF", "Fish")) {
ps <- prune_samples(sample_names(ps)[!grepl("GCSU21BBEB6", sample_names(ps))], ps)
} else {ps <- ps}
if (Datasetname %in% c("WF", "GCWF", "OEWF")) {
ps <- prune_samples(sample_names(ps)[!grepl("WFSU22MLEB", x =
sample_names(ps ))], ps)
} else {ps<-ps} #Missing physioligical data
################################
#Create Relative Abundance Data#
################################
ps_alpha_barplot <- ps %>%
tax_glom(taxrank = TaxLevel) %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>% # Transform to rel. abundance
psmelt() %>% # Melt to long format
filter(Abundance > 1) %>% # Filter out low abundance taxa
arrange(Genus) %>% # Sort data frame alphabetically by phylum
dplyr::arrange(desc(Abundance))
ps_alpha_barplot$ASV <- sub(".*:", "", ps_alpha_barplot$OTU)
ps_alpha_barplot$ASV <- sub('\\.', ' ', ps_alpha_barplot$ASV)
############################
#Create TotalAbundance Data#
############################
phylum_abundance <- ps_alpha_barplot %>%
dplyr::group_by(.data[[TaxLevel]]) %>%
dplyr::summarise(TotalAbundance = sum(Abundance))
ordered_levels <- phylum_abundance %>%
dplyr::arrange(desc(TotalAbundance)) %>%
pull(.data[[TaxLevel]])
ps_alpha_barplot$Taxa <- factor(ps_alpha_barplot[[TaxLevel]], levels = ordered_levels)
ps_alpha_barplot$SampleID <- factor(ps_alpha_barplot$SampleID, levels =
SSU_Samples[SSU_Samples %in% sample_names(pslist[[Dataset]])])
##########################
#Subset for Interest ASVs#
##########################
df <- ps_alpha_barplot #[ps_alpha_barplot$OTU %in% GroupOfInterest,]
##################################################
#Create Color Dataframe for Taxa from PhylaColors#
##################################################
taxonomy_df <- df[c("Phylum", "Taxa")]
#taxonomy_df$ASV <- sub(".*:", "", taxonomy_df$OTU)
taxonomy_df <- taxonomy_df[!duplicated(taxonomy_df$Taxa),]
#############################################
#Order to ordered_levels from Total Abudance#
order_index <- match(ordered_levels, taxonomy_df$Taxa)
taxonomy_df <- taxonomy_df[order_index, ]
print(head(taxonomy_df[c("Phylum", "Taxa")]))
taxonomy_df$Phylum2 <- ifelse(taxonomy_df$Phylum == "Proteobacteria",
paste0("Proteobacteria", ave(seq_along(taxonomy_df$Phylum),
taxonomy_df$Phylum, FUN = function(x) ifelse(seq_along(x) %% 2 == 0, 2, 1))),
taxonomy_df$Phylum)
unique_phyla <- unique(taxonomy_df$Phylum2)
# Assign base colors to unique phyla
#phylum_colors <- setNames(base_colors[1:length(unique_phyla)], unique_phyla)
# Print the mapping of phyla to colors
phylum_cols<- as.data.frame(phylum_colors)
phylum_cols$Phyla <- rownames(phylum_cols)
colored_taxonomy <- data.frame(
Taxa = taxonomy_df$Taxa,
Phylum2 = taxonomy_df$Phylum2,
Color = character(length(taxonomy_df$Taxa)), # Initialize an empty character vector
stringsAsFactors = FALSE)
# Generate shades and assign colors to each taxon
for (ii in seq_along(unique_phyla)) {
phylum <- unique_phyla[ii]
if (phylum %in% names(phylum_colors)) {
base_color <- phylum_colors[[phylum]]
} else {
base_color <- phylum_colors[["Other"]]
}
taxon_indices <- taxonomy_df$Phylum2 == phylum
num_variations <- sum(taxon_indices)
shades <- generate_shades(base_color, pmin(num_variations, 5))
#colored_taxonomy$Color[taxon_indices] <- shades
colored_taxonomy$Color[which(taxon_indices)[1:pmin(num_variations, 5)]]<- shades
}
colored_taxonomy$Color <- gsub("^$", NA, trimws(colored_taxonomy$Color))
colored_taxonomy$Color <- replace_na(colored_taxonomy$Color, "grey50")
##################
# #Some adjustments#
updateColor <- function(taxa, color) {
if (any(colored_taxonomy$Taxa == taxa)) {
colored_taxonomy[colored_taxonomy$Taxa == taxa,]$Color <- color
}
return(colored_taxonomy)
}
colored_taxonomy <- updateColor("Psychrobacter", "#FFCC00")
colored_taxonomy <- updateColor("Polynucleobacter", "#FF6600")
#updateColor("Photobacterium", "darkorange3")
# updateColor("Vibrio", "#FF9966")
#updateColor("Citrobacter", "#FFCC00")
#updateColor("Psychrobacter", "#FFCC00")
# updateColor("Enterobacter", "#993300")
#updateColor("Acinetobacter", "#FF9933")
Alpha_Diversity_df <- left_join(df, colored_taxonomy)
#################################
#Update to colored_taxonomy_Fish#
#################################
for (i in 1:nrow(Alpha_Diversity_df)) {
matching_row <- colored_taxonomy_Fish[colored_taxonomy_Fish$Taxa %in% Alpha_Diversity_df$Taxa[i], ]
if (nrow(matching_row) > 0) {
Alpha_Diversity_df$Color[i] <- matching_row$Color
}
}
color_mapping <- setNames(Alpha_Diversity_df$Color, Alpha_Diversity_df$Taxa)
levels(Alpha_Diversity_df$Taxa)[!levels(Alpha_Diversity_df$Taxa) %in% ordered_levels[1:34]] <-
"Other"
if (Datasetname == "WF") {
p <- ggplot(Alpha_Diversity_df,
aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
facet_grid(. ~ factor(Season, levels = SeasonOrder), drop = TRUE, scale = "free",
space = "free_x")
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in%
unique(sample_data(ps)$Season)]
} else {
p <- ggplot(Alpha_Diversity_df,
aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
facet_grid(. ~ factor(Reps, levels = RepOrder), drop = TRUE, scale = "free", space = "free_x")
COL <- col.Palette$col.Palette.Reps[names(col.Palette$col.Palette.Reps) %in%
unique(sample_data(ps)$Reps)]
}
p <- p +
atheme2 +
ylab("Relative Abundance [%] \n") + xlab("") +
labs(fill="") +
scale_fill_manual(values = color_mapping) +
#ggrepel::geom_text_repel(x = df$SampleID, y = df$Abundance, aes(label = Labels),
#inherit.aes = FALSE, min.segment.length = 0,nudge_y = 20,size=3, max.overlaps = Inf) +
#awhite +
theme(strip.text.y = element_text(angle = 0)) +
#theme(axis.text.x=element_blank(),
#axis.text.y=element_blank(),
#axis.title.y.left = element_blank(),
#axis.ticks.y =element_blank()
#) +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
#legend.background = element_rect(fill='transparent'), #transparent legend bg
#legend.box.background = element_rect(fill='transparent') #transparent legend panel
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=9),
legend.text=element_text(size=9),
axis.text.x.bottom = element_text(size= 9, face = "bold", angle = 45, hjust = 1),
strip.text.y.left = element_text(size = 9,face = "bold"),
axis.title.y.left = element_text(size=12,face = "bold"))
g <- ggplot_gtable(ggplot_build(p))
stripr <- which(grepl('strip-t', g$layout$name))
fills <- alpha(COL, 0.8)
k <- 1
for (iii in stripr) {
j <- which(grepl('rect', g$grobs[[iii]]$grobs[[1]]$childrenOrder))
g$grobs[[iii]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1}
Alpha_Diversity_BarPlot <- cowplot::plot_grid(g, labels = c(""), ncol = 1)
Alpha_Diversity_BarPlot <- plot_grid(Alpha_Diversity_BarPlot, ncol = 1, rel_heights = c(100))
ggsave(Alpha_Diversity_BarPlot, filename = FILENAME, path = pathPlots, device='png', dpi=300, width = WIDTH, height = 6)
plot(Alpha_Diversity_BarPlot )
}## Phylum Taxa
## 2 Actinobacteriota hgcI clade
## 3 Actinobacteriota CL500-29 marine group
## 5 Proteobacteria Limnohabitans
## 4 Verrucomicrobiota Persicirhabdus
## 55 Proteobacteria Polynucleobacter
## 28 Nitrospirota Nitrospira## Phylum Taxa
## 3 Bacteroidota Elizabethkingia
## 143 Proteobacteria Enterobacter
## 165 Proteobacteria Citrobacter
## 4 Proteobacteria Psychrobacter
## 35 Verrucomicrobiota Luteolibacter
## 133 Bacteroidota Chryseobacterium
## Phylum Taxa
## 8 Bacteroidota Elizabethkingia
## 78 Proteobacteria Enterobacter
## 77 Proteobacteria Citrobacter
## 4 Proteobacteria Polynucleobacter
## 95 Bacteroidota Chryseobacterium
## 99 Proteobacteria Psychrobacter
## Phylum Taxa
## 3 Bacteroidota Elizabethkingia
## 148 Proteobacteria Enterobacter
## 175 Proteobacteria Citrobacter
## 4 Proteobacteria Psychrobacter
## 35 Verrucomicrobiota Luteolibacter
## 135 Bacteroidota Chryseobacterium
## Phylum Taxa
## 8 Bacteroidota Elizabethkingia
## 79 Proteobacteria Enterobacter
## 76 Proteobacteria Citrobacter
## 4 Proteobacteria Polynucleobacter
## 99 Bacteroidota Chryseobacterium
## 102 Proteobacteria Psychrobacter
4.4 Core-Microbiome
Some thoughts from: McMurtrie, J., Alathari, S., Chaput, D. L., Bass, D., Ghambi, C., Nagoli, J., … & Tyler, C. R. (2021). Relationships between pond water and tilapia skin microbiomes in aquaculture ponds in Malawi. bioRxiv.McMurtrie, J., Alathari, S., Chaput, D. L., Bass, D., Ghambi, C., Nagoli, J., … & Tyler, C. R. (2021). Relationships between pond water and tilapia skin microbiomes in aquaculture ponds in Malawi. bioRxiv. Core microbiome analysis seeks to find taxa conserved between samples above a defined prevalence threshold. In this case we are interested in finding the core MB of the Fish Swab community which is present in fish across different Locations sites. Core calculations will be performed with the microbiome package. We will then visualise these core taxa as a heatmap using pheatmap. This will also reveal abundances across both Loc and fish samples, indicating if the fish core taxa are in high abundance in the Loc water and therefore potentially non-residents. Final visualisations for the finished figure requires some manual editing in Inkscape.
R Setup Libraries
Import phyloseq objects Weight up in choosing what taxa level to study. The full ASV set is large and requires a low core threshold as there is a reasonable degree of variation in ASVs e.g. likely different strains. However, in the core analysis we are most interested in the major types of organisms appearing to be common in either system. For instance we could have three individual ASVs of the same genus, each a key photosynthetic but only one is dominant per pond site, collectively though they are occupying the same niche. This is obviously a generalisation and not applicable to all scenarios. Most core studies I have seen group at genus level. I did explore the use of tip_gloming instead to use a phylogenetic informed grouping level, rather than arbitrary taxonomy, but this compromises resolution of taxonomic classifications that I am not willing to accept.
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE|ps_Fish", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")) {
ps <- pslist[[Dataset]]
Datasetname <- sub("ps_", "", Dataset)
pslist[[paste("Core", Datasetname, sep="_")]] <- list()
require(phyloseq)
core_taxa_standard <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
microbiome::core_members(., detection = 0.0001, prevalence = 90/100)
#ps_rel <- microbiome::transform(pslist[[Dataset]], "compositional")
#core_taxa_standard <- microbiome::core_members(ps_rel, detection = 0.0001, prevalence = 90/100)
#ps_core <- microbiome::core(ps_rel, detection = 0.0001, prevalence = .9)
ps_core <- prune_taxa(core_taxa_standard, ps)
core_taxa <- microbiome::taxa(ps_core)
tax_mat <- tax_table(ps_core)
tax_df <- as.data.frame(tax_mat)
tax_df$ASV <- rownames(tax_df)
# select taxonomy of only those OTUs that are core memebers based on the thresholds that were used.
core_taxa_class <- dplyr::filter(tax_df, rownames(tax_df) %in% core_taxa)
knitr::kable(head(core_taxa_class))
ps_rel_gen <- microbiome::aggregate_taxa(ps_core, "Genus")
# Check if any taxa with no genus classification. aggregate_taxa will merge all unclassified to Unknown
any(taxa_names(ps_rel_gen) == "Unknown")
# Remove Unknown
ps_rel_gen <- subset_taxa(ps_rel_gen, Genus!="Unknown")
library(RColorBrewer)
prevalences <- seq(.05, 1, .05)
detections <- round(10^seq(log10(1e-5), log10(.2), length = 10), 3)
# p1 <- microbiome::plot_core(ps_rel_gen,
# plot.type = "heatmap",
# colours = rev(brewer.pal(5, "RdBu")),
# prevalences = prevalences,
# detections = detections, min.prevalence = .8) +
# xlab("Detection Threshold (Relative Abundance (%))")
# p1 <- p1 + theme_bw() + ylab("ASVs")
# p1
pslist[[paste("Core", Datasetname, sep="_")]] <- core_taxa_class
}4.4.1 Plot Core
TaxLevel <- "ASV"
flipped <- T
for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])]) {
require(plyr); require(dplyr); require(ggrepel); require(cowplot); require(phyloseq)
Datasetname <- sub("ps_", "", Dataset)
ps <- pslist[[Dataset]]
ps2 <- ps #pslist[[paste(Dataset, "WF", sep="")]]
core_taxa_standard <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
microbiome::core_members(., detection = 0.0001, prevalence = 90/100)
GroupOfInterest <- core_taxa_standard
FILENAME <- paste(paste(save_name, "Core_BarPlot_noWF", Datasetname, sep="_"), ".png", sep="")
if (Datasetname %in% c("GC", "OE")) {
WIDTH <- 10 + length(sample_names(pslist[[Dataset]])) *0.12
} else if (Datasetname %in% c("WF")) {
WIDTH <- 5 + length(sample_names(pslist[[Dataset]])) *0.12
} else {
WIDTH <- 10 + length(sample_names(pslist[[Dataset]])) *0.12
}
################################
#Create Relative Abundance Data#
################################
ps_alpha_barplot <- ps2 %>%
#tax_glom(taxrank = TaxLevel) %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>% # Transform to rel. abundance
psmelt() %>% # Melt to long format
#filter(Abundance > 1) %>% # Filter out low abundance taxa
arrange(Genus) %>% # Sort data frame alphabetically by phylum
dplyr::arrange(desc(Abundance))
ps_alpha_barplot$ASV <- sub(".*:", "", ps_alpha_barplot$OTU)
#ps_alpha_barplot$ASV <- sub('\\.', ' ', ps_alpha_barplot$ASV)
############################
#Create TotalAbundance Data#
############################
# phylum_abundance <- ps_alpha_barplot %>%
# dplyr::group_by(.data[[TaxLevel]]) %>%
# dplyr::summarise(TotalAbundance = sum(Abundance))
# ordered_levels <- phylum_abundance %>%
# dplyr::arrange(desc(TotalAbundance)) %>%
# pull(.data[[TaxLevel]])
#
# ps_alpha_barplot$Taxa <- factor(ps_alpha_barplot[[TaxLevel]], levels = ordered_levels)
ps_alpha_barplot$SampleID <- factor(ps_alpha_barplot$SampleID, levels =
SSU_Samples[SSU_Samples %in% sample_names(ps2)])
SSU_Samples[SSU_Samples %in% sample_names(ps2)]
##########################
#Subset for Interest ASVs#
##########################
df <- ps_alpha_barplot #[ps_alpha_barplot$OTU %in% GroupOfInterest,]
#df <- df[df$OTU %in% GroupOfInterest,]
df$Abundance[!(df$OTU %in% GroupOfInterest)] <- 0
phylum_abundance <- df %>%
dplyr::group_by(.data[[TaxLevel]]) %>%
dplyr::summarise(TotalAbundance = sum(Abundance))
ordered_levels <- phylum_abundance %>%
dplyr::arrange(desc(TotalAbundance)) %>%
pull(.data[[TaxLevel]])
df$Taxa <- factor(df[[TaxLevel]], levels = ordered_levels)
Alpha_Diversity_df <- df
Alpha_Diversity_df$Colors <- NA
Alpha_Diversity_df$Reps <- factor(Alpha_Diversity_df$Reps, levels = RepOrder[RepOrder %in% Alpha_Diversity_df$Reps])
#################################
#Update to colored_taxonomy_Fish#
#################################
taxa_levels <- c("Genus", "Family", "Order", "Class", "Phylum")
taxa_level2 <- c("Taxa", "Phylum2")
Alpha_Diversity_df$Color <- "grey"
# Loop through each row of Alpha_Diversity_df
for (i in 1:nrow(Alpha_Diversity_df)) {
matching_color <- NA # Initialize matching color to NA
# Loop through each taxonomic level
for (level in taxa_levels) {
matching_row <- colored_taxonomy_Fish[colored_taxonomy_Fish$Taxa %in%
Alpha_Diversity_df[[level]][i], ]
# If there's a match, update matching_color and exit the loop
if (nrow(matching_row) > 0) {
matching_color <- matching_row$Color
break
}
}
# If there's no match at any level, try matching without numbers
if (is.na(matching_color)) {
matching_row <- colored_taxonomy_Fish[gsub("\\d", "", colored_taxonomy_Fish$Phylum2) %in%
Alpha_Diversity_df[["Phylum"]][i], ]
if (nrow(matching_row) > 0) {
matching_color <- matching_row$Color
}
}
# Assign the matching color to the corresponding row in Alpha_Diversity_df
Alpha_Diversity_df$Color[i] <- matching_color
}
color_mapping <- setNames(Alpha_Diversity_df$Color, Alpha_Diversity_df$Taxa)
color_mapping[["Candidatus Megaira"]] <- "purple"
color_mapping[["Acinetobacter.lwoffii"]] <- "#996600"
color_mapping[["Acinetobacter.johnsonii"]] <- "#663300"
color_mapping[["Shewanella"]] <- "#FF3366"
color_mapping[["Shewanella.baltica"]] <- "#FF0099"
color_mapping[["Shewanella.putrefaciens"]] <-"#FF0066"
color_mapping[["Photobacterium"]] <- "#FFF000"
color_mapping[["Aeromonas"]] <- "#FFCC00"
levels(Alpha_Diversity_df$Taxa)[!levels(Alpha_Diversity_df$Taxa) %in% head(unique(sub(".*:", "", GroupOfInterest)), 34) ] <-
"Other"
if (flipped == FALSE) {
if (Datasetname == "WF") {
p <- ggplot(Alpha_Diversity_df,
aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
facet_grid(. ~ factor(Season, levels = SeasonOrder), drop = TRUE, scale = "free",
space = "free_x")
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in%
unique(sample_data(ps2)$Season)]
} else {
p <- ggplot(Alpha_Diversity_df,
aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
facet_grid(. ~ factor(Reps, levels = RepOrder), drop = TRUE, scale = "free", space = "free_x")
COL <- col.Palette$col.Palette.Reps[names(col.Palette$col.Palette.Reps) %in% unique(sample_data(ps2)$Reps)]
}
p <- p +
atheme +
ylab("Relative Abundance [%] \n") + xlab("") +
labs(fill="") +
scale_fill_manual(values = color_mapping) +
ylim(0, 60) +
#ggrepel::geom_text_repel(x = df$SampleID, y = df$Abundance, aes(label = Labels),
#inherit.aes = FALSE, min.segment.length = 0,nudge_y = 20,size=3, max.overlaps = Inf) +
#awhite +
theme(strip.text.y = element_text(angle = 0)) +
#theme(axis.text.x=element_blank(),
#axis.text.y=element_blank(),
#axis.title.y.left = element_blank(),
#axis.ticks.y =element_blank()
#) +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
#legend.background = element_rect(fill='transparent'), #transparent legend bg
#legend.box.background = element_rect(fill='transparent') #transparent legend panel
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=9),
legend.text=element_text(size=9),
axis.text.x.bottom = element_text(size= 9, face = "bold", angle = 45, hjust = 1),
strip.text.y.left = element_text(size = 9,face = "bold"),
axis.title.y.left = element_text(size=12,face = "bold"))
g <- ggplot_gtable(ggplot_build(p))
stripr <- which(grepl('strip-t', g$layout$name))
fills <- alpha(COL, 0.8)
k <- 1
for (iii in stripr) {
j <- which(grepl('rect', g$grobs[[iii]]$grobs[[1]]$childrenOrder))
g$grobs[[iii]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1}
Alpha_Diversity_BarPlot <- cowplot::plot_grid(g, labels = c(""), ncol = 1)
Alpha_Diversity_BarPlot <- plot_grid(Alpha_Diversity_BarPlot, ncol = 1, rel_heights = c(100))
ggsave(Alpha_Diversity_BarPlot, filename = FILENAME, path = pathPlots, device='png', dpi=300,
width = WIDTH, height = 6)
} else if (flipped == TRUE) {
# New facet label names for supp variable
COL <- col.Palette$col.Palette.Reps[names(col.Palette$col.Palette.Reps) %in% unique(sample_data(ps2)$Reps)]
Short.labs <- gsub("[^0-9]", "", Alpha_Diversity_df$Reps)
names(Short.labs) <- Alpha_Diversity_df$Reps
p <- ggplot(Alpha_Diversity_df, aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
coord_flip() +
#facet_grid(rows = vars(Reps), drop = TRUE, scale = "free", space = "free_x", switch = "y")
facet_grid(rows = vars(Reps), drop = TRUE, scale = "free", space = "free_y", switch = "y",
labeller = labeller(Reps = Short.labs))
p <- p +
ylab("Relative Abundance [%] \n") + xlab("") +
labs(fill="") + #atheme+
scale_fill_manual(values = color_mapping) +
ylim(0, 80) +
theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank()) +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
) +
theme(axis.title.x.bottom = element_text(color="grey13", size=15, face ="bold"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=9),
legend.text=element_text(size=9),
axis.text.x.bottom = element_text(size= 15, face = "bold", angle = 45, hjust = 1),
#strip.text.y.left = element_text(size = 9,face = "bold"),
#axis.title.y.left = element_text(size=12,face = "bold")
) +
theme(legend.position = "none")
g <- ggplot_gtable(ggplot_build(p))
stripr <- which(grepl('strip-l', g$layout$name))
fills <- alpha(COL, 0.8)
k <- 1
for (iii in stripr) {
j <- which(grepl('rect', g$grobs[[iii]]$grobs[[1]]$childrenOrder))
g$grobs[[iii]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1}
Core_BarPlot <- cowplot::plot_grid(g, labels = c(""), ncol = 1)
Core_BarPlot <- plot_grid(Core_BarPlot, ncol = 1, rel_heights = c(100))
ggsave(Core_BarPlot, filename = FILENAME, path = pathPlots, device='png', dpi=300, width = 2.7, height = 10)
plot(Core_BarPlot)
}
}
#-
4.5 Bacterioplankton vs Fish (time and space)
Analyzing here the taxa that are shared between bacterioplankton and the fish gill mucus, in the first step we determine the “most shared” taxa by geometric mean of overlapping taxa between datasets In the second step we run mantel tests for the biological samples replicates and the bacterioplankton to see which groups show the strongest overall matrix correlation Finally we plot the shared taxa between bacterioplankton and fish mucus but excluding the fish gill core taxa which we saw before highly abundant in the fish but rare in the bacterioplankton
4.5.1 Most overlapping taxa
for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])]) {
Datasetname <- sub("ps_","", Dataset)
ps <- pslist[[paste("ps", Datasetname, sep="_")]]
SampleSums <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
phyloseq::otu_table() %>%
as.data.frame() %>%
rownames_to_column(var = "SampleID") %>%
dplyr::left_join(SAMDF16S, by = c("SampleID" = "SampleID")) %>%
dplyr::group_by(SampleID) %>%
dplyr::summarise(dplyr::across(rownames(tax_table(ps)), mean, na.rm = TRUE)) %>%
t() %>%
as.data.frame() %>%
`colnames<-`(.[1, ]) %>%
.[-1, ] %>%
stats::setNames(paste0("avg_", colnames(.))) %>%
mutate_all(as.numeric) %>%
round(., digits = 2) %>%
rownames_to_column(var="ASV")
SeasonSums <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
phyloseq::otu_table() %>%
as.data.frame() %>%
rownames_to_column(var = "SampleID") %>%
dplyr::left_join(SAMDF16S, by = c("SampleID" = "SampleID")) %>%
dplyr::group_by(Season) %>%
dplyr::summarise(dplyr::across(rownames(tax_table(ps)), mean, na.rm = TRUE)) %>%
t() %>%
as.data.frame() %>%
`colnames<-`(.[1, ]) %>%
.[-1, ] %>%
stats::setNames(paste0("avg_", colnames(.))) %>%
mutate_all(as.numeric) %>%
round(., digits = 2) %>%
rownames_to_column(var="ASV")
RepsSums <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
phyloseq::otu_table() %>%
as.data.frame() %>%
rownames_to_column(var = "SampleID") %>%
dplyr::left_join(SAMDF16S, by = c("SampleID" = "SampleID")) %>%
dplyr::group_by(Replicates) %>%
dplyr::summarise(dplyr::across(rownames(tax_table(ps)), mean, na.rm = TRUE)) %>%
t() %>%
as.data.frame() %>%
`colnames<-`(.[1, ]) %>%
.[-1, ] %>%
stats::setNames(paste0("avg_", colnames(.))) %>%
mutate_all(as.numeric) %>%
round(., digits = 2) %>%
rownames_to_column(var="ASV")
##########################
#Fish vs Bacterioplankton#
##########################
WF_ASVmeans <- pslist$ps_WF %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
phyloseq::otu_table() %>%
t() %>%
as.data.frame() %>%
dplyr::mutate(ASVmeans_WF = rowMeans(.)) %>%
dplyr::mutate(ASV = rownames(.)) %>%
dplyr::select(ASV, ASVmeans_WF)
Fish_ASVmeans <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
phyloseq::otu_table() %>%
t() %>%
as.data.frame() %>%
dplyr::mutate(ASVmeans_Fish = rowMeans(.)) %>%
dplyr::mutate(ASV = rownames(.)) %>%
dplyr::select(ASV, ASVmeans_Fish)
#Merge values on TAXON level
Fish_ASVmeans$Taxon <- sub("ASV\\d+:", "", Fish_ASVmeans$ASV)
Fish_Taxonmeans <- Fish_ASVmeans %>%
dplyr::group_by(Taxon) %>%
dplyr::summarize(Taxon_means_Fish = sum(ASVmeans_Fish)) %>%
dplyr::arrange(desc(Taxon_means_Fish))
WF_ASVmeans$Taxon <- sub("ASV\\d+:", "", WF_ASVmeans$ASV)
WF_Taxonmeans <- WF_ASVmeans %>%
dplyr::group_by(Taxon) %>%
dplyr::summarize(Taxon_means_WF = sum(ASVmeans_WF)) %>%
dplyr::arrange(desc(Taxon_means_WF))
#Combine Datasets
Overlapping_Fish <- Fish_Taxonmeans[Fish_Taxonmeans$Taxon %in% WF_Taxonmeans$Taxon,]
Overlapping_WF <- WF_Taxonmeans[WF_Taxonmeans$Taxon %in% Fish_Taxonmeans$Taxon,]
Overlapping <- Overlapping_Fish %>% left_join(Overlapping_WF)
###########################################################
#Take geometric mean of overlapping taxa in fish and Water#
print(paste("Most overlapping taxa between", Datasetname, "and Bakterioplankton"))
print(Overlapping %>%
dplyr::mutate(CombinedAbundance = sqrt(Taxon_means_Fish * Taxon_means_WF)) %>%
dplyr::arrange(desc(CombinedAbundance)))
#################################
#Overlapping ASVs betweem biomes#
#################################
#Combine Datasets
Overlapping_Fish_ASV <- Fish_ASVmeans[Fish_ASVmeans$ASV %in% WF_ASVmeans$ASV ,]
Overlapping_WF_ASV <- WF_ASVmeans[WF_ASVmeans$ASV %in% Fish_ASVmeans$ASV ,]
Overlapping_ASV <- Overlapping_Fish_ASV %>% left_join(Overlapping_WF_ASV)
# print(paste("Relative Abundance of Overlapping taxa between biomes in", Datasetname))
# print(sum(Overlapping$ASVmeans))
print(paste("Relative Abundance of Overlapping taxa between biomes in Bacterioplankton"))
print(sum(WF_ASVmeans[WF_ASVmeans$ASV %in% Overlapping_ASV$ASV,]$ASVmeans))
print(paste("Relative Abundance of Overlapping taxa between biomes in Bacterioplankton and", Datasetname))
RepsSums_BiomeSums <- RepsSums[RepsSums$ASV %in% Overlapping_ASV$ASV,] %>%
dplyr::summarize(across(starts_with("avg_"), sum, na.rm = TRUE))%>%
gather(., key = "Replicates", value = "SummedAbundance") %>%
dplyr::arrange(desc(SummedAbundance))
print(RepsSums_BiomeSums)
print(mean(RepsSums_BiomeSums$SummedAbundance))
RepsSums_BiomeSums_NoCore <- RepsSums[RepsSums$ASV %in% setdiff(Overlapping_ASV$ASV, pslist[[paste("Core", Datasetname, sep="_")]]$ASV),] %>%
dplyr::summarize(across(starts_with("avg_"), sum, na.rm = TRUE))%>%
gather(., key = "Replicates", value = "SummedAbundance") %>%
dplyr::arrange(desc(SummedAbundance))
print(RepsSums_BiomeSums_NoCore)
print(mean(RepsSums_BiomeSums_NoCore$SummedAbundance))
}## [1] "Most overlapping taxa between OE and Bakterioplankton"
## # A tibble: 247 × 4
## Taxon Taxon_means_Fish Taxon_means_WF CombinedAbundance
## <chr> <dbl> <dbl> <dbl>
## 1 Luteolibacter 4.09 1.58 2.54
## 2 Persicirhabdus 1.00 2.27 1.51
## 3 Flavobacterium 1.10 1.39 1.24
## 4 CL500-29 marine group 0.288 4.88 1.19
## 5 Elizabethkingia 16.7 0.0642 1.04
## 6 TRA3-20 0.549 1.44 0.889
## 7 Ilumatobacter 0.753 1.00 0.868
## 8 Halioglobus 0.468 1.46 0.828
## 9 Rhizobiales Incertae Sedis 0.400 0.986 0.628
## 10 Vicinamibacteraceae 0.262 1.45 0.617
## # ℹ 237 more rows
## [1] "Relative Abundance of Overlapping taxa between biomes in Bacterioplankton"
## [1] 48.71117
## [1] "Relative Abundance of Overlapping taxa between biomes in Bacterioplankton and OE"
## Replicates SummedAbundance
## 1 avg_OESP22MG 68.64
## 2 avg_OESU21MG 60.00
## 3 avg_OESP22BB 59.34
## 4 avg_OESP22ML 56.59
## 5 avg_OESU21BB 54.88
## 6 avg_OESP22SS 50.35
## 7 avg_OEAU21ML 49.69
## 8 avg_OESU21SS 49.61
## 9 avg_OESU21ML 49.01
## 10 avg_OESU21TW 47.51
## 11 avg_OEWI22SS 44.64
## 12 avg_OESP22TW 43.75
## 13 avg_OEWI22ML 43.43
## 14 avg_OEWI22MG 41.23
## 15 avg_OEWI22TW 40.13
## 16 avg_OEAU21TW 39.49
## 17 avg_OEWI22BB 36.43
## 18 avg_OEAU21SS 31.63
## 19 avg_OEAU21BB 28.21
## 20 avg_OEAU21MG 20.46
## [1] 45.751
## Replicates SummedAbundance
## 1 avg_OEAU21ML 47.47
## 2 avg_OESP22MG 40.21
## 3 avg_OEAU21TW 36.97
## 4 avg_OESU21ML 36.58
## 5 avg_OEAU21SS 29.07
## 6 avg_OESU21TW 27.65
## 7 avg_OESP22ML 26.26
## 8 avg_OEAU21BB 26.25
## 9 avg_OESU21MG 25.65
## 10 avg_OESU21SS 22.27
## 11 avg_OEWI22ML 20.85
## 12 avg_OESP22BB 20.57
## 13 avg_OESU21BB 19.60
## 14 avg_OEWI22MG 17.80
## 15 avg_OESP22SS 17.21
## 16 avg_OEWI22SS 16.63
## 17 avg_OEWI22TW 16.42
## 18 avg_OESP22TW 16.13
## 19 avg_OEAU21MG 15.28
## 20 avg_OEWI22BB 15.06
## [1] 24.6965
## [1] "Most overlapping taxa between GC and Bakterioplankton"
## # A tibble: 220 × 4
## Taxon Taxon_means_Fish Taxon_means_WF CombinedAbundance
## <chr> <dbl> <dbl> <dbl>
## 1 Polynucleobacter 4.40 0.825 1.91
## 2 Luteolibacter 1.09 1.58 1.31
## 3 Elizabethkingia 17.8 0.0642 1.07
## 4 Rhodobacteraceae 2.18 0.509 1.05
## 5 Persicirhabdus 0.438 2.27 0.997
## 6 Flavobacterium 0.582 1.39 0.901
## 7 TRA3-20 0.537 1.44 0.880
## 8 CL500-29 marine group 0.134 4.88 0.810
## 9 Ilumatobacter 0.531 1.00 0.729
## 10 Halioglobus 0.325 1.46 0.689
## # ℹ 210 more rows
## [1] "Relative Abundance of Overlapping taxa between biomes in Bacterioplankton"
## [1] 41.39703
## [1] "Relative Abundance of Overlapping taxa between biomes in Bacterioplankton and GC"
## Replicates SummedAbundance
## 1 avg_GCAU21TW 66.12
## 2 avg_GCAU21SS 58.72
## 3 avg_GCSP22BB 55.74
## 4 avg_GCAU21ML 54.60
## 5 avg_GCSU21BB 53.81
## 6 avg_GCSP22ML 50.78
## 7 avg_GCSP22TW 50.76
## 8 avg_GCWI22SS 50.74
## 9 avg_GCWI22BB 49.80
## 10 avg_GCSU21ML 46.99
## 11 avg_GCSP22SS 45.33
## 12 avg_GCAU21BB 42.48
## 13 avg_GCSU21TW 36.79
## 14 avg_GCWI22ML 35.65
## 15 avg_GCSU21SS 35.43
## [1] 48.916
## Replicates SummedAbundance
## 1 avg_GCSU21ML 32.26
## 2 avg_GCAU21TW 31.37
## 3 avg_GCWI22BB 24.46
## 4 avg_GCSU21TW 20.04
## 5 avg_GCWI22SS 17.25
## 6 avg_GCAU21ML 15.10
## 7 avg_GCSU21BB 13.81
## 8 avg_GCAU21SS 12.78
## 9 avg_GCSP22TW 10.66
## 10 avg_GCAU21BB 10.60
## 11 avg_GCSP22ML 10.31
## 12 avg_GCSP22BB 10.06
## 13 avg_GCSU21SS 9.63
## 14 avg_GCWI22ML 8.45
## 15 avg_GCSP22SS 8.24
## [1] 15.6684.5.2 Mantel tests with spatio-temporal resolution
####################################################
#Mantel test Bacterioplankton vs Fish per Replicate#
####################################################
Bacterioplankton_Mantel_list <- list()
for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])]) {
Datasetname <- sub("ps_","", Dataset)
Bacterioplankton_Mantel_list[[Datasetname]] <- list()
ps <- pslist[[paste("ps", Datasetname, sep="_")]]
RepsSums <- ps %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>%
phyloseq::otu_table() %>%
as.data.frame() %>%
rownames_to_column(var = "SampleID") %>%
dplyr::left_join(SAMDF16S, by = c("SampleID" = "SampleID")) %>%
dplyr::group_by(Replicates) %>%
dplyr::summarise(dplyr::across(rownames(tax_table(ps)), mean, na.rm = TRUE)) %>%
t() %>%
as.data.frame() %>%
`colnames<-`(.[1, ]) %>%
.[-1, ] %>%
stats::setNames(paste0("avg_", colnames(.))) %>%
mutate_all(as.numeric) %>%
round(., digits = 2) %>%
rownames_to_column(var="ASV")
RepsSums$Taxon <- sub("ASV\\d+:", "", RepsSums$ASV)
RepsSums_Taxonmeans <- RepsSums %>%
dplyr::group_by(Taxon) %>%
dplyr::summarize(across(starts_with("avg_"), sum, na.rm = TRUE))
RepsSums_Taxonmeans <- RepsSums_Taxonmeans %>%
dplyr::rowwise() %>%
dplyr::mutate(RowSum = sum(c_across(starts_with("avg_")), na.rm = TRUE)) %>%
dplyr::arrange(desc(RowSum))
merged_df <- inner_join(RepsSums_Taxonmeans, WF_Taxonmeans, by = "Taxon")
merged_df2 <- as.data.frame(merged_df)
rownames(merged_df2) <- merged_df2$Taxon
merged_df2 <- merged_df2[,grepl("avg_", colnames(merged_df2))]
colnames(merged_df2) <- sub("avg_", "", colnames(merged_df2))
TAXA <- as.data.frame(tax_table(ps))
TAXA$Taxa <- sub("ASV\\d+:", "", rownames(TAXA ))
TAXA <- TAXA[!duplicated(TAXA$Taxa),][c("Phylum", "Taxa")]
long_df<- merged_df2 %>%
pivot_longer(cols = everything(),
names_to = "Replicates",
values_to = "Abundance")
long_df <- cbind(Taxa = rownames(merged_df2), long_df)
long_df<- long_df %>%
left_join(SAMDF16S[!duplicated(SAMDF16S$Replicates),]) %>%
left_join(TAXA)
#long_df$Abundance[!(long_df$Taxa %in% GroupOfInterest)] <- 0
phylum_abundance <- long_df %>%
dplyr::group_by(.data[["Taxa"]]) %>%
dplyr::summarise(TotalAbundance = sum(Abundance))
ordered_levels <- phylum_abundance %>%
dplyr::arrange(desc(TotalAbundance)) %>%
pull(.data[["Taxa"]])
long_df $Taxa <- factor(long_df[["Taxa"]], levels = ordered_levels)
calculate_correlations <- function(df, columns, comparison_column) {
correlations <- sapply(columns, function(col) {
stats::cor(df[[col]], df[[comparison_column]], use = "complete.obs",
method = c("pearson"))
})
return(correlations)
}
# Get the columns you want to compare
columns_to_compare <- names(RepsSums_Taxonmeans)[grep("^avg_", names(RepsSums_Taxonmeans))]
# Calculate the correlations
correlations <- calculate_correlations(merged_df, columns_to_compare, "Taxon_means_WF")
correlations %>%
as.data.frame() %>%
dplyr::arrange(desc(correlations))
#############################
#Mantel Test on every column#
#############################
mantel_results <- lapply(columns_to_compare, function(COLUM) {
# Create distance matrices for the column and Taxon_means_WF
dist_matrix_col <- as.matrix(dist(merged_df[[COLUM]], method = "euclidean"))
dist_matrix_wf <- as.matrix(dist(merged_df$Taxon_means_WF, method = "euclidean"))
# Perform the Mantel test
mantel_test <- mantel(dist_matrix_col, dist_matrix_wf, permutations = 999)
return(c(mantel_test$statistic, mantel_test$signif))
})
mantel_results_df <- as.data.frame(do.call(rbind, mantel_results))
colnames(mantel_results_df) <- c("Mantel_Statistic", "P_Value")
mantel_results_df <- cbind(Column = columns_to_compare, mantel_results_df)
# Order the results by Mantel statistic
mantel_results_df <- mantel_results_df %>%
dplyr::arrange(desc(Mantel_Statistic))
mantel_results_df$Replicates <- sub("avg_", "", mantel_results_df$Column)
mantel_results_df <- mantel_results_df %>%
dplyr::left_join(SAMDF16S[SAMDF16S$Species %in% Datasetname & !duplicated(SAMDF16S$Replicates),])
mantel_results_df$Reps <- factor(mantel_results_df$Reps, levels = RepOrder)
Bacterioplankton_Mantel_list[[Datasetname]] <- mantel_results_df
}
Bacterioplankton_Mantel <- rbind(Bacterioplankton_Mantel_list$OE, Bacterioplankton_Mantel_list$GC)
#####################
#Plot Mantel Results#
#####################
Bacterioplankton_Mantel$Season <- factor(Bacterioplankton_Mantel$Season, levels=SeasonOrder)
RepOrder2 <-
c("SU 633", "SU 651", "SU 665", "SU 692", "SU 713",
"AU 633", "AU 651", "AU 665", "AU 692", "AU 713",
"WI 633", "WI 651", "WI 665", "WI 692", "WI 713",
"SP 633", "SP 651", "SP 665", "SP 692", "SP 713")
Bacterioplankton_Mantel$Reps <- factor(Bacterioplankton_Mantel$Reps, levels=RepOrder2)
Bacterioplankton_Mantel$Species <- factor(Bacterioplankton_Mantel$Species, levels =c("OE", "GC"))
Mantel_plot <- Bacterioplankton_Mantel %>%
ggplot(., aes(x=Species, y=Reps, fill=Mantel_Statistic )) +
geom_tile() +
scale_fill_gradient2(
low = "#FFFFFF",
high = "#006000",
limit = c(-0.2,1)) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, size = 15, face ="bold"),
#axis.title.x.bottom = element_text(size=15,face = "bold"),
#axis.title.y.left = element_text(size=15,face = "bold"),
strip.text.x = element_text(size = 15, face = "bold"),
strip.text.x.bottom = element_text(size = 15,face = "bold"),
strip.text.y.left = element_text(size = 15,face = "bold"),
strip.text.y.right = element_text(size = 15,face = "bold"),
axis.title.x.bottom = element_blank(),
axis.title.y.left = element_blank(),
axis.text.x.bottom = element_text(face = "bold", angle = 45, hjust = 1),
axis.text.y.left = element_text(face = "bold", size = 15),
axis.text.y.right = element_text(face = "bold", size = 15),
legend.title = element_text(size = 15, face ="bold"),
legend.text = element_text(size = 15,face = "bold"),
plot.title = element_text(size = 15, face = "bold"),
plot.subtitle = element_text(size = 15),
plot.caption = element_text(size = 15)) +
labs(#title = "Correlation relationship between avg.dist Bray-Curtis matrices and varibles", y = "Sample kind",
fill="corr") +
labs(fill = "Mantel's r") +
#theme(legend.box.background = element_rect(color = "black"))
theme(legend.position="right", legend.key.width = unit(1, "cm"))
Mantel_plot <- Mantel_plot + geom_text(vjust=0.76,hjust=0.49, aes(label = paste0(c("","*")[(abs(P_Value) <= .05)+1],
c("","*")[(abs(P_Value) <= .001)+1], c("","*")[(abs(P_Value) <= .001)+1])))
# Mantel_plot <- Mantel_plot + facet_grid(. ~ factor(Season, levels = SeasonOrder), drop = TRUE, scale = #"free", space = "free_x")
Mantel_plot <- Mantel_plot + facet_grid(vars(Season), scale="free")
Mantel_plot <- cowplot::plot_grid(Mantel_plot, labels = c(""), ncol = 1)
Mantel_plot <- plot_grid(Mantel_plot, ncol = 1, rel_heights = c(0.02, 0.05, 0.98))
plot(Mantel_plot)
ggsave(Mantel_plot, filename = paste(paste("Mantel_test_Bacterioplankton_FishReps",
sep="_"),".png", sep=""), path = pathPlots , device='png', dpi=300, width = 4,height = 7)4.5.3 Stacked barplot of overlapping taxa without core
colored_taxonomy_Fish <- pslist$Optics$colored_taxonomy_Fish
#################
#GroupofInterest#
#################
GroupOfInterest <- setdiff(Overlapping_ASV$ASV, pslist$Core_Fish$ASV)
Comparison_Name <- "Biomes_Overlapping_Taxa_noCore"
TaxLevel <- "ASV"
flipped <- T
for (Dataset in names(pslist)[grepl("ps_GCWF|ps_OEWF", names(pslist))]) {
require(plyr); require(dplyr); require(ggrepel); require(cowplot); require(phyloseq)
Datasetname <- sub("ps_", "", Dataset)
ps <- pslist[[Dataset]]
GroupOfInterest <- GroupOfInterest
FILENAME <- paste(paste(save_name, "Alpha_BarPlot", Comparison_Name, Datasetname, sep="_"), ".png", sep="")
if (Datasetname %in% c("GC", "OE")) {
WIDTH <- 10 + length(sample_names(pslist[[Dataset]])) *0.12
HEIGHT <- 10 #length(sample_names(pslist[[Dataset]])) *0.08
} else if (Datasetname %in% c("WF")) {
WIDTH <- 5 + length(sample_names(pslist[[Dataset]])) *0.12
HEIGHT <- 10 #length(sample_names(pslist[[Dataset]])) *0.08
} else {
WIDTH <- 10 + length(sample_names(pslist[[Dataset]])) *0.12
HEIGHT <- 10 #length(sample_names(pslist[[Dataset]])) *0.08
}
################################
#Create Relative Abundance Data#
################################
ps_alpha_barplot <- ps %>%
#tax_glom(taxrank = TaxLevel) %>%
transform_sample_counts(function(x) {x/sum(x)*100}) %>% # Transform to rel. abundance
psmelt() %>% # Melt to long format
#filter(Abundance > 1) %>% # Filter out low abundance taxa
arrange(Genus) %>% # Sort data frame alphabetically by phylum
dplyr::arrange(desc(Abundance))
ps_alpha_barplot$ASV <- sub(".*:", "", ps_alpha_barplot$OTU)
#ps_alpha_barplot$ASV <- sub('\\.', ' ', ps_alpha_barplot$ASV)
############################
#Create TotalAbundance Data#
############################
# phylum_abundance <- ps_alpha_barplot %>%
# dplyr::group_by(.data[[TaxLevel]]) %>%
# dplyr::summarise(TotalAbundance = sum(Abundance))
# ordered_levels <- phylum_abundance %>%
# dplyr::arrange(desc(TotalAbundance)) %>%
# pull(.data[[TaxLevel]])
#
# ps_alpha_barplot$Taxa <- factor(ps_alpha_barplot[[TaxLevel]], levels = ordered_levels)
ps_alpha_barplot$SampleID <- factor(ps_alpha_barplot$SampleID, levels =
SSU_Samples[SSU_Samples %in% sample_names(pslist[[Dataset]])])
SSU_Samples[SSU_Samples %in% sample_names(pslist[[Dataset]])]
##########################
#Subset for Interest ASVs#
##########################
df <- ps_alpha_barplot #[ps_alpha_barplot$OTU %in% GroupOfInterest,]
#df <- df[df$OTU %in% GroupOfInterest,]
df$Abundance[!(df$OTU %in% GroupOfInterest)] <- 0
phylum_abundance <- df %>%
dplyr::group_by(.data[[TaxLevel]]) %>%
dplyr::summarise(TotalAbundance = sum(Abundance))
ordered_levels <- phylum_abundance %>%
dplyr::arrange(desc(TotalAbundance)) %>%
pull(.data[[TaxLevel]])
df$Taxa <- factor(df[[TaxLevel]], levels = ordered_levels)
Alpha_Diversity_df <- df
Alpha_Diversity_df$Colors <- NA
Alpha_Diversity_df$Reps <- factor(Alpha_Diversity_df$Reps, levels = RepOrder[RepOrder %in%
Alpha_Diversity_df$Reps])
#################################
#Update to colored_taxonomy_Fish#
#################################
taxa_levels <- c("Genus", "Family", "Order", "Class", "Phylum")
taxa_level2 <- c("Taxa", "Phylum2")
Alpha_Diversity_df$Color <- "grey"
# Loop through each row of Alpha_Diversity_df
for (i in 1:nrow(Alpha_Diversity_df)) {
matching_color <- NA # Initialize matching color to NA
# Loop through each taxonomic level
for (level in taxa_levels) {
matching_row <- colored_taxonomy_Fish[colored_taxonomy_Fish$Taxa %in%
Alpha_Diversity_df[[level]][i], ]
# If there's a match, update matching_color and exit the loop
if (nrow(matching_row) > 0) {
matching_color <- matching_row$Color
break
}
}
# If there's no match at any level, try matching without numbers
if (is.na(matching_color)) {
matching_row <- colored_taxonomy_Fish[gsub("\\d", "", colored_taxonomy_Fish$Phylum2) %in%
Alpha_Diversity_df[["Phylum"]][i], ]
if (nrow(matching_row) > 0) {
matching_color <- matching_row$Color
}
}
# Assign the matching color to the corresponding row in Alpha_Diversity_df
Alpha_Diversity_df$Color[i] <- matching_color
}
color_mapping <- setNames(Alpha_Diversity_df$Color, Alpha_Diversity_df$Taxa)
color_mapping[["Candidatus Megaira"]] <- "purple"
color_mapping[["Acinetobacter.lwoffii"]] <- "#996600"
color_mapping[["Acinetobacter.johnsonii"]] <- "#663300"
color_mapping[["Shewanella"]] <- "#FF3366"
color_mapping[["Shewanella.baltica"]] <- "#FF0099"
color_mapping[["Shewanella.putrefaciens"]] <-"#FF0066"
color_mapping[["Photobacterium"]] <- "#FFF000"
color_mapping[["Aeromonas"]] <- "#FFCC00"
levels(Alpha_Diversity_df$Taxa)[!levels(Alpha_Diversity_df$Taxa) %in% head(unique(sub(".*:", "", GroupOfInterest)), 34) ] <-
"Other"
if (flipped == FALSE) {
if (Datasetname == "WF") {
p <- ggplot(Alpha_Diversity_df,
aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
facet_grid(. ~ factor(Season, levels = SeasonOrder), drop = TRUE, scale = "free",
space = "free_x")
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in%
unique(sample_data(ps)$Season)]
} else {
p <- ggplot(Alpha_Diversity_df,
aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
facet_grid(. ~ factor(Reps, levels = RepOrder), drop = TRUE, scale = "free", space = "free_y")
COL <- col.Palette$col.Palette.Reps[names(col.Palette$col.Palette.Reps) %in% unique(sample_data(ps)$Reps)]
}
p <- p +
atheme +
ylab("Relative Abundance [%] \n") + xlab("") +
labs(fill="") +
scale_fill_manual(values = color_mapping) +
ylim(0, 100) +
#ggrepel::geom_text_repel(x = df$SampleID, y = df$Abundance, aes(label = Labels),
#inherit.aes = FALSE, min.segment.length = 0,nudge_y = 20,size=3, max.overlaps = Inf) +
#awhite +
theme(strip.text.y = element_text(angle = 0)) +
#theme(axis.text.x=element_blank(),
#axis.text.y=element_blank(),
#axis.title.y.left = element_blank(),
#axis.ticks.y =element_blank()
#) +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
#legend.background = element_rect(fill='transparent'), #transparent legend bg
#legend.box.background = element_rect(fill='transparent') #transparent legend panel
) +
theme(axis.title.x.bottom = element_text(color="grey13"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=9),
legend.text=element_text(size=9),
axis.text.x.bottom = element_text(size= 9, face = "bold", angle = 45, hjust = 1),
strip.text.y.left = element_text(size = 9,face = "bold"),
axis.title.y.left = element_text(size=12,face = "bold"))
g <- ggplot_gtable(ggplot_build(p))
stripr <- which(grepl('strip-t', g$layout$name))
fills <- alpha(COL, 0.8)
k <- 1
for (iii in stripr) {
j <- which(grepl('rect', g$grobs[[iii]]$grobs[[1]]$childrenOrder))
g$grobs[[iii]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1}
Alpha_Diversity_BarPlot <- cowplot::plot_grid(g, labels = c(""), ncol = 1)
Alpha_Diversity_BarPlot <- plot_grid(Alpha_Diversity_BarPlot, ncol = 1, rel_heights = c(100))
ggsave(Alpha_Diversity_BarPlot, filename = FILENAME, path = pathPlots, device='png', dpi=300,
width = WIDTH, height = 6)
} else if (flipped == TRUE) {
# New facet label names for supp variable
COL <- col.Palette$col.Palette.Reps[names(col.Palette$col.Palette.Reps) %in% unique(sample_data(ps)$Reps)]
Short.labs <- gsub("[^0-9]", "", Alpha_Diversity_df$Reps)
names(Short.labs) <- Alpha_Diversity_df$Reps
p <- ggplot(Alpha_Diversity_df, aes(x = SampleID, y = Abundance, fill = factor(Taxa))) +
geom_bar(stat = "identity") +
coord_flip() +
#facet_grid(rows = vars(Reps), drop = TRUE, scale = "free", space = "free_x", switch = "y")
facet_grid(rows = vars(Reps), drop = TRUE, scale = "free", space = "free_x", switch = "y",
labeller = labeller(Reps = Short.labs))
p <- p +
ylab("Relative Abundance [%] \n") + xlab("") +
labs(fill="") + #atheme+
scale_fill_manual(values = color_mapping) +
ylim(0, 100) +
theme(axis.title.y=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank()) +
theme(
panel.background = element_rect(fill='transparent'), #transparent panel bg
plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
) +
theme(axis.title.x.bottom = element_text(color="grey13", size=15, face ="bold"),
strip.text = element_text(color = "black", face= "bold")) +
theme(
legend.title=element_text(size=9),
legend.text=element_text(size=9),
axis.text.x.bottom = element_text(size= 15, face = "bold", angle = 45, hjust = 1),
#strip.text.y.left = element_text(size = 9,face = "bold"),
#axis.title.y.left = element_text(size=12,face = "bold")
) +
theme(legend.position = "none")
g <- ggplot_gtable(ggplot_build(p))
stripr <- which(grepl('strip-l', g$layout$name))
fills <- alpha(COL, 0.8)
k <- 1
for (iii in stripr) {
j <- which(grepl('rect', g$grobs[[iii]]$grobs[[1]]$childrenOrder))
g$grobs[[iii]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1}
Alpha_Diversity_BarPlot <- cowplot::plot_grid(g, labels = c(""), ncol = 1)
Alpha_Diversity_BarPlot <- plot_grid(Alpha_Diversity_BarPlot, ncol = 1, rel_heights = c(100))
ggsave(Alpha_Diversity_BarPlot, filename = FILENAME, path = pathPlots, device='png', dpi=300, width = 2.7, height = HEIGHT)
plot(Alpha_Diversity_BarPlot)
}
}
-
5 Alpha diversity
first steps of alpha diversity analysis the dataset shows overall little patterns and high multicollinearity between influencing variables, therefore we decided to skip modeling alpha diversity
5.1 Metrics & Stats
###################
#Phylum Statistics#
###################
Phylum_Statistics_list <- list()
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")) {
individual <- pslist[[Dataset]] %>% tax_glom(taxrank = "Phylum") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
#dplyr::filter(Abundance > 0.0005) %>%
dplyr::arrange("Phylum")
individual_abund <- as.data.frame(individual %>%
dplyr::select(Sample, SampleID, Replicates2, Phylum, Abundance) %>%
dplyr::group_by(Phylum, Sample) %>%
dplyr::summarise(TotalAbundance = (sum(Abundance))*100)) %>%
dplyr::ungroup()%>%
dplyr::left_join(individual %>%
dplyr::select(Sample, SampleID, Replicates2, Phylum, Abundance),
by = c("Sample", "Phylum"))
individual_abund_SD <- individual_abund %>%
dplyr::group_by(Phylum) %>%
dplyr::summarise(SD = round(sd(Abundance)*100, digits = 1))
merged_species <- merge_samples(pslist[[Dataset]], "Species")
# Use psmelt to obtain a long-format data.frame
merged_species <- merged_species %>% tax_glom(taxrank = "Phylum") %>%
transform_sample_counts(function(x) {x/sum(x)}) %>%
psmelt() %>%
#dplyr::filter(Abundance > 0.0005) %>%
dplyr::arrange("Phylum")
merged_species_abund <- as.data.frame(merged_species %>%
dplyr::select(Sample, SampleID, Phylum, Abundance) %>%
dplyr::group_by(Phylum, Sample) %>%
dplyr::summarise(TotalAbundance = round(sum(Abundance)*100,digits = 1)))
print(Dataset)
data <- merged_species_abund %>%
left_join(individual_abund_SD) %>%
dplyr::arrange(desc(TotalAbundance))
print(data)
Phylum_Statistics_list[[Dataset]] <- data
}## [1] "ps_OE"
## Phylum Sample TotalAbundance SD
## 1 Proteobacteria OE 51.0 10.3
## 2 Bacteroidota OE 27.7 11.9
## 3 Actinobacteriota OE 6.8 5.3
## 4 Verrucomicrobiota OE 5.3 7.6
## 5 Firmicutes OE 4.0 7.1
## 6 Deinococcota OE 2.4 2.3
## 7 Chloroflexi OE 0.7 1.2
## 8 Acidobacteriota OE 0.4 0.5
## 9 Patescibacteria OE 0.4 0.4
## 10 Desulfobacterota OE 0.3 0.6
## 11 Cyanobacteria OE 0.2 0.4
## 12 Dependentiae OE 0.2 0.9
## 13 Nitrospirota OE 0.2 0.4
## 14 Bdellovibrionota OE 0.1 0.2
## 15 Myxococcota OE 0.1 0.2
## 16 Planctomycetota OE 0.1 0.2
## 17 Armatimonadota OE 0.0 0.1
## 18 Campylobacterota OE 0.0 0.3
## 19 Fusobacteriota OE 0.0 0.3
## 20 Gemmatimonadota OE 0.0 0.0
## 21 Spirochaetota OE 0.0 0.1
## [1] "ps_GC"
## Phylum Sample TotalAbundance SD
## 1 Proteobacteria GC 57.1 13.2
## 2 Bacteroidota GC 25.0 9.1
## 3 Actinobacteriota GC 5.9 4.8
## 4 Firmicutes GC 5.0 13.3
## 5 Verrucomicrobiota GC 2.6 2.3
## 6 Deinococcota GC 2.2 2.2
## 7 Acidobacteriota GC 0.5 0.7
## 8 Chloroflexi GC 0.5 0.8
## 9 Cyanobacteria GC 0.2 0.3
## 10 Desulfobacterota GC 0.2 0.3
## 11 Nitrospirota GC 0.2 0.3
## 12 Patescibacteria GC 0.2 0.3
## 13 Planctomycetota GC 0.2 0.2
## 14 Armatimonadota GC 0.1 0.1
## 15 Bdellovibrionota GC 0.0 0.1
## 16 Campylobacterota GC 0.0 0.1
## 17 Dependentiae GC 0.0 0.1
## 18 Fusobacteriota GC 0.0 0.0
## 19 Gemmatimonadota GC 0.0 0.0
## 20 Myxococcota GC 0.0 0.1
## 21 Spirochaetota GC 0.0 0.1
## [1] "ps_WF"
## Phylum Sample TotalAbundance SD
## 1 Proteobacteria WF 41.1 5.4
## 2 Bacteroidota WF 19.2 7.4
## 3 Actinobacteriota WF 19.1 8.8
## 4 Verrucomicrobiota WF 8.4 3.3
## 5 Acidobacteriota WF 2.3 1.8
## 6 Nitrospirota WF 1.9 1.4
## 7 Desulfobacterota WF 1.5 1.1
## 8 Cyanobacteria WF 1.4 1.9
## 9 Gemmatimonadota WF 1.4 0.6
## 10 Chloroflexi WF 1.2 0.9
## 11 Myxococcota WF 0.7 0.4
## 12 Planctomycetota WF 0.4 0.6
## 13 Bdellovibrionota WF 0.3 0.2
## 14 Firmicutes WF 0.3 0.4
## 15 Campylobacterota WF 0.2 0.4
## 16 Armatimonadota WF 0.1 0.1
## 17 Deinococcota WF 0.1 0.1
## 18 Dependentiae WF 0.1 0.0
## 19 Patescibacteria WF 0.1 0.0
## 20 Spirochaetota WF 0.1 0.1
## 21 Caldisericota WF 0.0 0.0
## 22 Calditrichota WF 0.0 0.0
## 23 Cloacimonadota WF 0.0 0.0
## 24 Deferrisomatota WF 0.0 0.0
## 25 Hydrogenedentes WF 0.0 0.0
## 26 Latescibacterota WF 0.0 0.0
## 27 Methylomirabilota WF 0.0 0.0write.csv2(do.call(rbind, Phylum_Statistics_list), file.path(path_Output_16S,paste(paste(save_name,"Phylum_Statistics_list",Date, sep="_"), ".csv", sep="")))5.1.1 P:B Ratio
Ratios of Proteobacteria and Bacteroidetes
Rosado, D., Pérez-Losada, M., Severino, R., Cable, J., & Xavier, R. (2019). Characterization of the skin and gill microbiomes of the farmed seabass (Dicentrarchus labrax) and seabream (Sparus aurata). Aquaculture, 500, 57-64.
5.2 Alpha Indices
https://f1000research.com/articles/5-1492
alphadiv_list <- list()
for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE|ps_Fish|ps_SSU", names(pslist)[!grepl("WF", names(pslist))])]) {
require(phyloseq)
require(vegan)
require(plyr)
require(dplyr)
require(ggrepel)
require(cowplot)
Datasetname <- sub("ps_", "", Dataset)
a <- length(alphadiv_list)
Samples <- sample_names(pslist[[Dataset]])
ps <- merge_phyloseq(pslist$ps_OE, pslist$ps_GC, pslist$ps_WF)
ps <- prune_samples(sample_names(ps) %in% Samples, ps)
rowSums(otu_table(ps))
min(rowSums(otu_table(ps)))
max(rowSums(otu_table(ps)))
mean(rowSums(otu_table(ps)))
#Richness vs SequencingDepth#
ggplot(data = data.frame("total_reads" = phyloseq::sample_sums(ps),
"observed" = phyloseq::estimate_richness(ps, measures = "Observed")[, 1]),
aes(x = total_reads, y = observed)) +
geom_point() +
geom_smooth(method="lm", se = FALSE) +
labs(x = "\nTotal Reads", y = "Observed Richness\n")
if (Datasetname %in% c("GC", "OE", "GCWF", "OEWF", "Fish")) {
print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
print(which(rowSums(otu_table(ps)) < 5000))
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
} else if (Datasetname %in% c("SSU")) {
print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
print(which(rowSums(otu_table(ps)) < 5000))
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
} else if (Datasetname %in% c("WF")) {
print(paste("No Waterfiler Samples removed ;", Dataset, sep = ""))
print(rowSums(otu_table(ps)))
ps_Filter <- ps
}
ps_Filter_rare <- phyloseq::rarefy_even_depth(ps_Filter , rngseed = 123, replace = FALSE)
alphadiv_list[[a+1]] <- ps_Filter_rare
names(alphadiv_list)[[a+1]] <- paste("ps_rare", Datasetname, sep="_")
vegan::rarecurve(as.data.frame(otu_table(ps_Filter)), step=100, lwd=2, ylab="ASVs", label=F)
abline(v=(min(rowSums(as.data.frame(otu_table(ps_Filter))))))
vegan::rarecurve(as.data.frame(otu_table(ps_Filter_rare)), step=100, lwd=2, ylab="ASVs", label=F)
adiv <- data.frame(
"Observed" = phyloseq::estimate_richness(ps_Filter_rare, measures = "Observed"),
"Chao1" = phyloseq::estimate_richness(ps_Filter_rare, measures = "Chao1"),
"Shannon" = phyloseq::estimate_richness(ps_Filter_rare, measures = "Shannon"),
"InvSimpson" = phyloseq::estimate_richness(ps_Filter_rare, measures = "InvSimpson"),
"LOC" = phyloseq::sample_data(ps_Filter_rare)$LOC,
"Season" = phyloseq::sample_data(ps_Filter_rare)$Season,
"SampleID" = phyloseq::sample_data(ps_Filter_rare)$SampleID
)
alphadiv_list[[a+2]] <- adiv
names(alphadiv_list)[[a+2]] <- paste("adiv", Datasetname, sep="_")
}## [1] "Samples removed for low frequency below 5000 seqs in;ps_SSU"
## GCAU21TWFL2 GCAU21MLEB1 WFSU21SSFL WFWI22MGEB
## 183 184 229 232
## [1] "Samples removed for low frequency below 5000 seqs in;ps_Fish"
## GCAU21TWFL2 GCAU21MLEB1
## 183 184
## [1] "Samples removed for low frequency below 5000 seqs in;ps_OE"
## named integer(0)
## [1] "Samples removed for low frequency below 5000 seqs in;ps_GC"
## GCAU21TWFL2 GCAU21MLEB1
## 45 46
pslist$AlphaDiversity <- alphadiv_list
alphadiv_list <- pslist$AlphaDiversity5.3 Data Exploration
DataExploration_list <- list()
for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])]) {
require(phyloseq)
require(vegan)
require(plyr)
require(dplyr)
require(ggrepel)
require(cowplot)
require(ggordiplots)
Datasetname <- sub("ps_", "", Dataset)
ps <- pslist[[Dataset]]
DataExploration_list[[paste(Datasetname)]] <- list()
GroupingVariables <- c("Season", "LOC")
Traits <- c(
"Salinity","O2_catch","SecchiDepth", "Temperature","pH_catch","Age",
"FCF", "HSI", "SSI", "Length", "Weigth", "GSI", "FillLevel",
"NH4", "NO2", "NO3", "O2", "PO4", "pH", "SPM",
GroupingVariables)
print(paste("Samples removed for low frequency below 5000 seqs in;", Datasetname, sep = ""))
print(which(rowSums(otu_table(ps)) < 5000))
ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
'%!in%' <- function(x,y)!('%in%'(x,y))
Include <- Traits[Traits %!in% GroupingVariables]
df_Traits <- data.frame(sample_data(ps_Filter))
df_Traits <- df_Traits[names(df_Traits) %in% Traits]
df_Include <- df_Traits[names(df_Traits) %in% Include]
#Just to inspect how well FGG-Data fit to your measured values while fishing:
#df_Include[c("SPM", "SecchiDepth")]
#df_Include[c("pH_catch", "pH")]
#df_Include[c("O2_catch", "O2")]
print("NAs in TraitData")
print(table(is.na(df_Include)))
print((df_Include[!complete.cases(df_Include), ]))
df_Include <- df_Include[complete.cases(df_Include), ]
df_Traits <- df_Traits[complete.cases(df_Traits), ]
dim(df_Traits)
#df_Include <- vegan::decostand(df_Include, method='standardize')
dim(df_Include)
ps_Filter <- prune_samples(sample_names(ps_Filter)[sample_names(ps_Filter) %in%
rownames(df_Include[complete.cases(df_Include), ])], ps_Filter)
adiv <- pslist$AlphaDiversity[[paste("adiv", Datasetname, sep="_")]]
adiv <- adiv[rownames(adiv) %in% rownames(df_Include),]
model_df <- left_join(df_Include %>% mutate(SampleID = rownames(.)), adiv %>% mutate(SampleID = rownames(.)))
model_df$fSeason <- factor(model_df$Season, level = SeasonOrder)
model_df$fLOC <- factor(model_df$LOC, level = LOCOrder)
model_df$Species <- substr(model_df$SampleID, 1, nchar(model_df$SampleID) - 9)
model_df$Replicates <- substr(model_df$SampleID, 1, nchar(model_df$SampleID) - 3)
model_df$fReplicates <- factor(model_df$Replicates, levels = VariableOrder[VariableOrder %in% model_df$Replicates])
rownames(model_df) <- model_df$SampleID
print(Datasetname)
head(model_df)
DataExploration_list[[paste(Datasetname)]][[paste("ps_Filter")]] <- ps_Filter
DataExploration_list[[paste(Datasetname)]][[paste("model_df")]] <- model_df
}## [1] "Samples removed for low frequency below 5000 seqs in;OE"
## named integer(0)
## [1] "NAs in TraitData"
##
## FALSE
## 2622
## [1] Temperature pH_catch Salinity O2_catch SecchiDepth Length
## [7] FillLevel Age HSI SSI GSI FCF
## [13] NH4 NO2 NO3 O2 PO4 pH
## [19] SPM
## <0 rows> (or 0-length row.names)
## [1] "OE"
## [1] "Samples removed for low frequency below 5000 seqs in;GC"
## GCAU21TWFL2 GCAU21MLEB1
## 45 46
## [1] "NAs in TraitData"
##
## FALSE
## 1634
## [1] Temperature pH_catch Salinity O2_catch SecchiDepth Length
## [7] FillLevel Age HSI SSI GSI FCF
## [13] NH4 NO2 NO3 O2 PO4 pH
## [19] SPM
## <0 rows> (or 0-length row.names)
## [1] "GC"pslist$DataExploration_list <- DataExploration_list 5.3.1 Violin Plots
5.3.1.1 Alpha Between species
5.3.1.2 AlphaDiv
###################
#Variable Violin#
###################
alphadiv_violin <- list()
for (Dataset in names(alphadiv_list)[grepl("adiv", names(alphadiv_list))]) {
require(plyr)
require(ggrepel)
require(cowplot)
require(EnvStats)
Variables <- c("Observed","Chao1.Chao1","Shannon", "InvSimpson")
Datasetname <- sub("adiv_", "", Dataset)
adiv <- alphadiv_list[[Dataset]]
for (ii in Variables) {
alphadiv_violin_length <- length(alphadiv_violin)
FILENAME <- paste(paste(save_name, Type, Datasetname, ii, "Violin", sep="_"),".png", sep="")
ggg <- paste(ii)
adiv$LOC <- factor(adiv$LOC, levels=LOCOrder)
p <- ggplot(adiv, aes(x=LOC,y=adiv[[ii]], colour=LOC)) +
#geom_boxplot(aes(fill=LOC)) + atheme +
geom_violin(aes(fill=LOC)) +
geom_jitter(aes(color = LOC)) +
scale_color_manual(values= ggplot2::alpha(col.Palette$col.Palette.LOC, 1)) +
scale_fill_manual(values= ggplot2::alpha(col.Palette$col.Palette.LOC, 0.4)) +
#geom_text_repel(aes(label = SampleID), data = adiv, size=2.5) +
facet_grid(. ~ factor(Season, levels=SeasonOrder), drop=TRUE,scale="free",space="free_x") +
labs(x = "", y = noquote(paste(ggg, " [", noquote(AlphaIndices[ii]), "]", sep = ""))) + #"")
#ggtitle(paste(gg, ggg, "over all sampling sites and seasons (Summer, Winter, Spring,)"))
theme(legend.position = "none") +
guides(fill=guide_legend(title="Species")) +
#stat_n_text(size = 4, color = "grey20")+ #white
theme(
panel.background = element_rect(fill='transparent'),
plot.background = element_rect(fill='transparent', color=NA),
#panel.grid.major = element_blank(),
#panel.grid.minor = element_blank(),
legend.background = element_rect(fill='transparent'),
legend.box.background = element_rect(fill='transparent') ) +
theme(
axis.title.y.left = element_text(size=10,face = "bold"),
strip.text.y.left = element_text(size = 10,face = "bold"),
axis.text.y.left = element_text(face = "bold"),
legend.title = element_text( size = 6),
legend.text = element_text(size = 6,face = "bold"),
plot.title = element_text(size = 10, face = "bold"),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 9),
axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
strip.text.x = element_text(face = "bold"))
COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in% adiv$Season]
g <- ggplot_gtable(ggplot_build(p))
stripr <- which(grepl('strip-t', g$layout$name))
fills <- alpha(COL, 0.4)
k <- 1
for (iii in stripr) {
j <- which(grepl('rect', g$grobs[[iii]]$grobs[[1]]$childrenOrder))
g$grobs[[iii]]$grobs[[1]]$children[[j]]$gp$fill <- fills[k]
k <- k+1}
prow <- cowplot::plot_grid(g, labels = c(""), ncol = 1)
title <- ggdraw() + draw_label_themeRK(paste(Datasetname, ggg, "per Location"), element = "plot.title",x = 0.05, hjust = 0, vjust = 0.5) #draw_label_themeRKwhite
alphadiv_violin_plot<- cowplot::plot_grid(prow, ncol = 1, rel_heights = c(0.05, 0.98)) #title,
plot(alphadiv_violin_plot)
ggsave(alphadiv_violin_plot, filename = FILENAME, path = pathPlots, device='png', dpi=300, width = 5,
height = 2.8)
alphadiv_violin[[alphadiv_violin_length+1]] <- alphadiv_violin_plot
names(alphadiv_violin)[[alphadiv_violin_length+1]] <- paste(Datasetname, ii, sep="_")
}
}
##############
#Summary Plot#
##############
cowplot::plot_grid(alphadiv_violin$OE_Observed, alphadiv_violin$GC_Observed,
labels = c("A", "B"), ncol = 1) -> part_1
cowplot::plot_grid(alphadiv_violin$OE_Shannon, alphadiv_violin$GC_Shannon,
labels = c("C", "D"), ncol = 1) -> part_2
cowplot::plot_grid(part_1, part_2, labels = c("", ""), rel_widths = c(1,1), ncol = 2) -> part_3
ggsave(part_3, filename = paste(paste(Species, Year, Season, "AlphaDiv_Violin", Date, sep="_"),".png", sep=""), path = pathPlots , device='png', dpi=300, width = 13,height = 6)
part_3
5.3.1.3 Physiology
5.3.1.3 Physio line plot
#-
Save Data
###########
#Save Data#
###########
saveRDS(pslist, file.path(path_Output_16S, paste(paste(save_name, "ps_16S_merge_Filter_ASV_besttax", Date, sep="_"), ".rds", sep="")))